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Chapter 1 
INTRODUCTION 


1.1 Structural Dynamics 

Today, every airplane, booster, and spacecraft In production 
has been subjected to a very thorough dynamic analysis. However, it Is 
not only the aerospace Industry that is interested in structural 
dynamics. The Interest has spread through a range of Industries which 
manufacture everything from automobiles to typewriters, and construc- 
tion companies which build everything from bridges to high-rise office 
buildings. 

The reason for wanting a dynamic analysis may vary from company 
to company, but certain requirements of the analysis remain constant. 
It must be accurate. It must be low cost, and It must be completed 
quickly. Today's structural dynamlclsts are addressing each of these 
requirements. 

One of the most popular approaches to satisfying these require- 
ments Involves the use of substructures. 

1.2 Reasons for Substructuring 

Structural systems have become so complex today that no single 
engineer, and often no single company, is capable of efficiently 
designing and analyzing the entire structure alone. The structures are 
too large and too complex to handle on the system level, and are 
therefore broken down Into smaller structures (l.e. substructures). 
The three primary reasons for substructuring are: 
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1. Different groups or companies specialize In one particular part 0 * 
a structure, and are the most qualified to design that part. 
Therefore the structure Is separated into substructures, and each 
substructure Is given to the group which specializes in its design. 

2. When a single group or company is given responsibility for the 
entire structure, the project is often more easily managed if the 
structure is separated into substructures. Each substructure may 
then be designed and verified independently before assembling all 
the substructures to verify the superstructure. 

3. From a computational point of view, it is often the case that a 
coupled substructure analysis will yield a more efficient analysis 
of the structure than will a single analysis of the entire struc- 
ture. 


Having discussed the reasons for substructuring, the actual 
dynamic analysis will now be considered. General analysis techniques 
will be discussed first, and substructure analysis techniques will be 
discussed subsequently. 

1.3 Current Response Analyses 

The dynamic analysis of a structure may include several analy- 
ses including a modal analysis, a shock spectrum analysis, or a re- 
sponse analysis to a particular dynamic load. The primary focus of 
this paper will be on the latter. 
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Currently, there are two basic types of analyses used to obtain 
the system response to dynamic load. They may be classified as 
(1) time-domain analyses, and (2) frequency-domain analyses. 

Time-domain analysis techniques commonly attempt to uncouple 
the system equations and then Integrate them directly. This uncoupling 
of the equations of motion is accomplished by using a modal transforma- 
tion. Since the transformation matrix must be orthogonal with respect 
to the mass, stiffness, and damping matrices, it is often found by 
solving the eigenvalue problem which arises from the free vibration 
problem and subsequently assuming a form of damping which leaves the 
modes uncoupled. 

Frequency-domain analysis techniques use either Laplace or 
Fourier transformations to transform the differential eouations of 
motion to the frequency domain, where they become algebraic equations. 
The system of linear algebraic equations may then be solved for the 
system response. This frequency-domain response may then be inverse 
transformed to obtain the response in the time domain. 

1.4 Substructure Coupling in the Frequency Domain 

When utilizing a frequency-domain procedure for coupling 
substructures, only the dynamic stiffness matrices for each substruc- 
ture are required. The system equations are then obtained by using a 
direct stiffness approach to couple the substructure dynamic stiffness 
matrices. This set of linear equations may then be solved for the 
system response. 
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There are some instances where the frequency-domain techniques 
are more efficient than time-domain techniques. However, there are 
some instances where this is reversed. Since each of these methods is 
capable of giving either exact or approximate solutions, this paper 
will frequently compare the efficiency of the two procedures, and 
consequently will often discuss the cost of a particular operation. 
This cost is most accurately defined as the amount of standardized 
computer time required for that particular operation. 

1.5 Current Frequency Domain Analysis 

The best-known frequency-domain analysis is the harmonic or 
steady-state analysis. This frequency-domain technique has been used 
for many years to determine the steady-state response of structures 
subjected to periodic loads, and has been implemented in several of the 
large general purpose finite element packages (Ref. 1, 2). An equally 
well-known frequency-domain analysis is the determination of the 
spectral response of a structure to random loads (Ref. 3, 4). 

Little work has been done in the area of determining the 
transient response of structures using frequency-domain techniques. It 
is surmised that this is because of the problems encountered when 
transforming problems between the discrete time and frequency domains, 
and because of the expense of frequency-domain analysis of a large 
unreduced system. Each of these topics will be addressed in subsequent 
chapters of this thesis. 

An area of frequency-domain analysts which is becoming popular 
among analysts is the representation of substructures by frequency- 
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dependent matrices. These methods are often denoted as impedance 
methods, and various papers have been published on the topic by Chopra, 
Poelart, Payne, Geering, and Will lams and Wittrlck (Ref. 3, 5, 6, 7, 
8 ). 

Chopra uses frequency-domain methods to perform earthquake 
analyses of concrete dams. He uses frequency-domain representations of 
the soil and the dam which are derived from finite element models to 
obtain his system equations of motion. He reduces the number of 
degrees of freedom through the use of Rltz vectors. 

Poelart presents a distributed-element approach to the fre- 
quency-domain method, where an element is represented by its exact 
Impedance matrix. His paper also presents a discussion of the problems 
associated with using finite elements In a dynamic analysis, as he 
compares the distributed element .approach to the finite element ap- 
proach. 

The paper by Payne presents an approach by which two finite 
element substructure models may be coupled using the Impedance matrix 
of each substructure. He uses modal expansions of each substructure to 
determine the respective Impedance matrices. His paper gives parti- 
cular emphasis to problems with statical ly-determinate interfaces. 
Payne also provides for a reduction in the computational effort when 
structural modifications are performed on one of the substructures. 

The work by Geering Is primarily concerned with the coupling of 
substructures represented by frequency-dependent matrices, and the 
reduction in the number of degrees of freedom required for an analysis. 
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Each of these topics is discussed in subsequent chapters, as much of 
this paper is based on the work of Geering. 

The computation of natural frequencies using frequency-domain 
techniques has been Investigated by Williams and Wittrick. In their 
work, each of the elements of a structure is represented by its dynamic 
stiffness matrix. Their procedure for determining natural frequencies 
Is discussed in subsequent chapters also. 

Using these papers to establish the current state of the art In 
frequency-domain analysis, the research goals for this thesis were then 
determined. 

1.6 Research Goals 

The primary goal of this research project was to obtain a clear 
understanding of the strengths and weaknesses of substructure coupling 
in the frequency domain. These qualities were then to be evaluated and 
used to determine where the method is particularly applicable, and 
where it Is not. 

The secondary goal was to determine the generality of the 
method. This was to be accomplished by determining If a substructured 
model which had been created in the frequency domain could be used to 
obtain the system natural frequencies. 

It Is not the purpose of the research project to demonstrate 
new analytical techniques, although some have been presented. The 
purpose of the paper is to provide an analyst with sufficient informa- 
tion about the strengths and weaknesses of substructure coupling in the 
frequency domain. The analyst should then be able to make a decision 
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regarding the usefulness of the method for a particular application. 
An attempt has also been made to provide sufficient detail to permit 
implementation of the Geering method if so desired. 

1.7 Organization of the Thesis 

This thesis has been organized into eight chapters. The first 
chapter provides justification and background for the work presented in 
subsequent chapters. 

Chapter 2 describes the system equations of motion, and the 
analytical tools necessary to transform the equations to the frequency 
domain. The chapter discusses both steady state and transient analy- 
sis. Chapter 3 proceeds with a discussion of the implementation of 
frequency-domain analysis on a digital computer, and Chapter 4 dis- 
cusses Geering's method of coupling substructures in the frequency 
domain (Ref. 3). Chapters 2, 3 and 4 are included as support for the 
primary research goal. 

The support for the secondary research goal is included in 
Chapters 5 and 6, which outline procedures which may be used to deter- 
mine the natural frequencies of undamped and damped structures respec- 
tively. Chapter 5 is based on work by Williams and Wittrlck (Ref. 6) 
In their search for natural frequencies of continuous systems, while 
Chapter 6 presents some original work on determining the natural 
frequencies of damped systems. 

Chapter 7 presents examples to support the equations of Chap- 
ters 2 through 6, and Chapter 8 is the final chapter which includes 
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some general conments on substructure coupling in the frequency domain. 
Directions for further research are also presented in Chapter 8. 



Chapter 2 

FREQUENCY-DOMAIN ANALYSIS 


2.1 Equations of Motion 

Consider a finite element model of a structure with n active 
degrees of freedom.* The model may be formed using either lumped mass 
matrices or consistent mass matrices. The damping may be proportional , 
nonproportional, hysteretic, or any other linear form. In other words, 
the analyst is free to model the structure in various ways, as long as 
the equations of motion for the system may be written as 

[m] (x( t ) } + [c] (x(t)} + [k] (x(t)} ■ (f(t)} (2.1) 

[m], [c], and Tk] are time-invariant matrices commonly denoted 
as the mass, damping, and stiffness matrices respectively. The 
(x(t)} , (x ( t) ) , and (x(t)> vectors represent the acceleration, 
velocity, and displacement of each of the n degrees of freedom. 
Likewise, the (f (t) } vector represents the load applied at each of 
the n degrees of freedom. The initial displacement and initial 
velocity vectors will be written as 


(x(t=0)> = 

<v 

(2.2a) 

{x(t=0)> = 

<v 

(2.2b) 


^ Equations throughout this thesis will include only active degrees of 
freedom, which are those degrees of freedom which have not been 
constrained to have zero displacement. 
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Equation (2.1) represents a system of n linear second order 
differential equations. However, by applying carefully selected 
transformations, the equations may be transformed into a system of n 
linear algebraic equations. 

This chapter will examine three such transform pairs which 
allow the equations to be transformed to the new domain, and the 
solution to be Inverse transformed back to the time domain. Each 
transformation will be examined to determine (1) whether it gives a 
steady state or transient solution, (2) whether there are restrictions 
on the type of excitation, and (3) whether the system is permitted to 
be damped or undamped. 

2.2 Fourier Integral Transforms 

The Fourier Integral transform Is one of the most popular 
transformations because of the ease with which it may be discretized 
and Implemented on a digital computer. The unilateral Fourier Integral 
transform will be used in this paper, since all excitation and re- 
sponses will be assumed to be Identically zero before some instant In 
time (t=0 for this paper). 

Given a general function y(t) , the unilateral Fourier 
Integral transform pair may be written 


Y(«) « f 

Jo 

y(t) e" jwt dt 

(2.3a) 

ti 

4-> 

f* Y(w)e Jait du 

J «QO 

(2.3b) 


where co Is a real variable, often denoted as the frequency variable. 
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The equation of motion In the frequency domain may be obtained 
by applying the transformation given by Eq. (2.3a) to Eq. (2.1). 
Appendix A.l gives the details of the transformation which yields 


[-u 2 [m] + ju [c] + [k]] {X((u)} 

= {F(u)> + [m] {x Q > + [jo, [m] + [c]] {x Q } 


(2.4) 


or In shorthand notation, 


[G(u)] (X(o>) > a {F(d))) (2.5) 

where the dynamic stiffness matrix, [G(g>)] , Is given by 

[G(u)] = -u 2 [m] + joj [c] + [k] (2.6) 


and 

[F(u)] = {F(u)J + [m] {x 0 } + [ju> [ml + [c]]{x Q ) (2.7) 

where (F(a») > is the unilateral Fourier Integral transform of the 

forcing function {f (t) > , and {x Q > and {x Q } are the Initial 

displacement and velocity vectors respectively. 

The response spectrum may then be obtained by solving the 
system of linear equations In Eq. (2.5). Thus, 

(X(w) } = [H(o>)] {F(a>)} (2.8) 

where the dynamic flexibility matrix, [H(w)] , is given by 

[HU)] = CSU)]* 1 


(2.9) 
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The time-domain response may then be obtained by inverse transforming 
{X(w)} using Eq. (2.3b). 

Careful consideration must be given to the restrictions and 
assumptions made in Appendix A.l, as they determine the class of 
problems which may be solved using the Fourier integral transform 
approach. In short, they may be stated: 

- The forcing function (f ( t ) > must be of 
finite duration, and therefore nonperiodic. 

- The system must be damped. 

If these conditions are satisfied, the Inverse transform of the 
spectrum given by Eq. (2.8) will yield the exact, or as It Is often 
called, the transient solution to Eq. (2.1). 

In summary, the unilateral Fourier integral transform approach 
may be used to solve for the transient response of damped systems 
subjected to finite duration excitation. 

2.3 Fourier Integral Transforms with Convergence Functions 

An investigation of the Fourier transformabllity conditions in 
Appendix A.l reveals that the class of problems which may be worked 
using the Fourier integral transform approach is limited because of the 
condition that the Integral, 

f |y(t)| dt (2.10) 

j 0 

must converge. This convergence condition necessitated the restriction 
of finite-duration excitation and response. 
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Therefore, In order that the Fourier integral transform ap- 
proach may be used for systems not satisfying the convergence condi- 
tion, Eq. (2.1) will be multiplied by a yet unknown function, b(t) , 

b(t) [fm] {x> + [c] {x} + [k] (x)] = b(t) (f ( t) > (2.11) 

where b(t) Is any function which will cause Eq. (2.11) to satisfy the 
convergence conditions when transformed. The unilateral Fourier 
Integral transform may now be applied to Eq. (2.11) and written 


b(t) [[m] 


(X} + [c] {x> + 


[k] (x)] 


-jut 


dt 



b(t) {f(t)> e‘ jut dt 


( 2 . 12 ) 


Now assume that the convergence function, b(t) , is given by 

b(t) = e" at (2.13) 

where the convergence factor, a , is any value which will cause 
Eq. (2.11) to satisfy the convergence conditions when transformed. 

Equation (2.13) may now be substituted Into equation (2.12) to 

yield 


” [[m] {5} + fc] (x) + [k] <x>] e‘ at e' Jot dt = 
Jn J 



--at -jut 

e e 


(2.14) 


dt 
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An examination of the transformabll ity of Eq. (2.14) reveals that the 
integrals 


f [O] (J) + [c] (x) «• [k] <x)] e' at 


dt 


and 




dt 


(2.15) 


(2.16) 


must converge. It has been shown that for forcing functions of 
exponential order a Q , (i.e. for those functions which have the 
property that there Is a real number, a Q , such that 

(f(t)> e“ at = 0 when a > a Q (2.17) 

and with the limit not existing when a < a Q *») Eq. (2.16) will 
converge for a > o Q . (Ref. 9). It has also been shown that since 
the response of stable linear systems Is at least of exponential order 
0, Eq. (2.15) will always converge for a > 0 , and will converge for 
a s 0 when the system and excitation is of the type described In 
Section 2.2. 

In light of these guarantees for convergence, the restrictions 
which required the excitation and response to be of finite duration may 
be lifted. It Is Interesting to note, however, that the selection of a 
'large* convergence factor will cause the integrands of Eq. (2.12) to 
approach zero 'very fast,' simulating a finite duration excitation and 


response. 



Equation (2.14) may now be simplified and written 
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f [w {x} ♦ [c] {X} + [k] {x}] e" st dt 


f <f(t)} 
J 0 


e“ st dt 


(2.18) 


where the complex frequency variable, s , is given by 


s * a + j&> 


(2.19) 


At this point, Eq. (2.18) may be integrated by parts as de- 
tailed in Appendix A. 2, and written 


[w s' + Hi + Tk]] (X(s)) = 
CF(s)> + Cb]{x # > * [[hi] s ♦ W]<x 0 ) 


( 2 . 20 ) 


or in shorthand notation as 


[G(s)3 (X(s)} * (F(s)} 


( 2 . 21 ) 


where 


CG(s)] - W $* + [cl s + [kj 


( 2 . 22 ) 


and 


(F(s)} * (F(s)l + [m] {x_> + [[ml s + Fcj] (x } (2.23) 


It has been shown in Appendix A. 2 that 
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F(s) 


= f f(t) e“ st 
JO 


dt 


(2. 24) 


or 


• r 

Jn 


-at -jut 


F(s) = f(t) e e" J '"'' dt 
f O 


(2.25) 


which is simply the unilateral Fourier integral transform of (f ( t ) } , 


where 


(f(t)> = (f(t)> e 


■at 


(2.26) 


As in the previous section, Eq. (2.1) may be solved for 
(X(s)} , yielding 


where 


(X(s)> = [H(s)1 { F(s ) > 


[H(s)] = rG(s )]' 1 


(2.27) 


(2.28) 


Referencing Appendix A. 2 once again, reveals that 


{X(s ) } = f {x} e” st 


dt 


(2.29) 


or 


{X(s ) } 


f(x> e at 
J 0 


5 jut dt (2.30) 


which is simply the unilateral Fourier integral transform of (x) , 


where 
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{x} = {x} e" at (2.31) 

Therefore, the Inverse unilateral Fourier Integral transfom of {X(s ) > 
does not yield {x> , but instead yields {x> . However, (x > may be 
recovered using the relationship 

{x> = {x} e at (2.32) 

At this point, it should be clear that when the convergence 
function, e" at , is used in conjunction with the unilateral Fourier 
Integral transformation, the transform is equivalent to the Laplace 
transformation, where the complex frequency variable, s , Is simply 
the Laplace variable. The detailed discussion of the convergence 
function was undertaken In this section, (1) to emphasize why the 
Laplace transform removes some of the restrictions of the Fourier 
transform, (2) to emphasize that unilateral Fourier Integral transfor- 
mations may be used to perform forward and Inverse Laplace transforms 
once a convergence factor has been selected, and (3) to give some 
physical Insight into what exactly the convergence function does. In 
order to aid in the selection of the convergence factor. 

Therefore, the unilateral Fourier Integral transformation, used 
in conjunction with a convergence function, e” at , may be used to 
solve for the transient response of a damped or undamped system sub- 
jected to any excitation of exponential order. 
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2.4 Complex Fourier Series 

In the previous sections, methods of obtaining the transient 
response of a system were discussed. However, there are instances when 
only the steady-state response to a periodic excitation is desired. 

For this class of problems, the equations of motion may be 
transformed to the frequency domain using the complex Fourier series 
transform pair. 


where 


y(t) • 



Y H> 


jw. t 
e * 


(2.33) 


Y <"k> 


f T+Tl 

y(t) e 

T 




dt 


(2.34) 


•k 



(2.35) 


and Tj Is the fundamental period of y(t) . 

The details of the transformation are presented In Appendix 
A. 3. The frequency domain equations of motion are 


[-coi* [m] + ju. [c] + [kl ] {X(w|J} = 


{F(a> k )} 


(2.36) 


k = 


where (F(u k )} is the forcing spectrum determined by a complex Fourier 
analysis of (f ( t) } . As in the previous sections, Eq. (2.36) may be 
simplified and solved for (XU^)} 
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[G(« k )] 

rS 

3 

U_ 

II 

rS 

3 

X 

k = 

(2.37) 

3 

X 

= CH(- k )l {F(« k )> 

k = . . ,® 

(2.38) 


where 

[H(oj| ( )] = [G(to k )] -1 k = — ,...,» (2.39) 

The time-domain response may then be obtained using the Inverse trans- 
form relation given in Eq. (2.33). 

An examination of the transformation detailed in Appendix A. 3 
will reveal that the time-domain solution obtained is only a particu- 
lar, or steady-state, solution to the equations of motion given by Eq. 
( 2 . 1 ). 

Therefore, Chapter 2 has presented two transformations which 
may be used to obtain the transient response to a limited class of 
problems. The chapter has also presented a well-known method of 
obtaining the steady-state response of a system. Now, Chapter 3 will 
Investigate the computational considerations in the implementation of 
these procedures. 



Chapter 3 

COMPUTATIONAL CONSIDERATIONS 

3.1 Discretizing the Problem 

In Chapter 2, three different sets of frequency domain equa- 
tions of motion were presented. Those presented in Eqs. (2.4) and 
(2.20) were functions of the continuous variables o> and s respec- 
tively, and the equations of motion presented In Eq. (2.36) were 
functions of the discrete variable . It Is the purpose of this 
section to consolidate the three sets of equations into one set of 
discretized equations to be used In the discussion of the frequency- 
domain computational model. 

The first step in the consolidation Is to recognize that Eq. 
(2.4) Is identical to Eq. (2.20) when the convergence factor is identi- 
cally zero. Therefore, since Eq. (2.4) is contained within Eq. (2.20), 
the latter equation will be used henceforth. 

It has been shown by several authors, that the Fourier integral 
transform and the inverse Fourier integral transform, such as those 
required in the evaluation of Eq. (2.20), may be approximated by the 
discrete Fourier transform (DFT) pair (Ref. 4, 10, 11). 

Therefore, the discretization of Eq. (2.20) will begin by using 
the DFT to evaluate {F(s)> at k discrete frequencies. The 
transformation and the resulting force spectrum may be written as 

((F k ) ; k=0,...,K-l) = DFT ({f(t)>* e" at ) (3.1) 


20 
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ip 

where { > is equivalent to (F( ) > . The { > notation indicates 
that only discrete values of {f ( t ) > are used as required by the DR 
algorithm. 

Since the force spectrum, as given by Eq. (3.1), is only 
defined at K discrete frequencies, Eq. (2.20) may only be evaluated 
at those K frequencies. Thus 

[O] s£ + [c] s k + [k] ] (X k > 

(3.2) 

* {F k > + [m] {x Q } + [[m] s k + [c]] {x Q } 
k=0,...,K-l 

Equation (3.2) represents K sets of simultaneous equations 
which must be solved in order to obtain the discrete response spectrum 
({X k > ; k=0,...,K-l) . 

This discrete response spectrum may then be inverse transformed 

to obtain the discrete function, {x(t)>* e" at , and Eq. (2.31) may be 

* 

used to determine (x(t)} . Thus, 

(x(t)>* = e at [IDR ((X k > ; k=0,...,K-l)] (3.3) 

★ 

where the l ) notation again represents a discrete time function. 

Having discretized Eq. (2.20), the consolidation of Eqs. (3.2) 
and (2.36) into a single set of equations will continue. The next step 
In the consolidation Is to recognize the fact that functions of the 
variable . a) k may be represented as functions of s k with a 
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convergence factor of zero. Equation (2.36) may therefore be written 
as 


[[m] s£ + [c] s k + [kl] (X k ) = {F k > k— (3.4) 

Since it has been shown that the DFT may also be used to 
approximate the complex Fourier series (Ref. 4), the transformation of 
the periodic forcing function f f (t) > may be written as 

(fF k > ; k a 0,...,K-l) = DFT ({f(t)}*) (3.5) 

which is equivalent to Eq. (3.1) when the convergence factor is zero. 

Again, since the force spectrum is only defined at K discrete 
frequencies, Eq. (3.4) may only be evaluated at those frequencies. 
Thus, 


[[m] s£ + [c] s k + Tk]] {X k > = {F k > k=0,. . .K-l (3.6) 

Now, Eqs. (3.2) and (3.6) may both be written in shorthand 
notation as 



[6 k ] (X k > = {F k > k=0,...,K-l 

(3.7) 

where 



tv ■ 

[m] s| + [c] s k + [k] k=0,...,K-l 

(3.8) 

tF k> ■ ( V + 

[m] (x Q > + [fm] s k + [c]] (x 0 > k=0,...,K-l 

(3.9) 


and 
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({F k > ; k=0». . . ,K-1) = DFT ({f(t)>* e" at ) (3.10) 

The frequency-domain response spectrum may be found using the relation- 
ships 


(X k > = [H k ] {F k > k=0,...,K-l (3.11) 

where 

[H k ] = [G^' 1 k=0,. . . ,K-1 (3.12) 

and finally, the time-domain response may be found at discrete values 
of t from 


{x(t) }* = e at [IDFT ({X k } ; k=0,. . . ,K-1)] (3.13) 

Subject to the following modifications, Eqs. (3.7) through 
(3.13) represent the single computational model which may be imple- 
mented to obtain the response of the systems described in Sections 2.2, 
2.3, and 2.4. 

- To determine the transient response of damped systems to 
finite duration excitation 

* Set a = 0 

- To determine the transient response of either damped or 
undamped systems 

’ Select a proper value for the convergence factor (Ref. 


Section 3.3) 
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- To determine the steady-state response of either damped or 
undamped systems to a period input 

* Set a = 0 

* Set {x Q } = {x Q } = {0} (This Is simply a method of 
getting the algorithm to Ignore the initial conditions, 
since they are not used in steady-state analysis). 

Therefore, the solution procedures have been consolidated into 
a single set of equations which is easily adapted for each desired 
response analysis. In an attempt to ease the burden of notation in the 
remainder of this paper, Eqs. (3.7) and (3.11) will be written as 


[G] (X) = (F> 

(3.14) 

(X) = [H] {F} 

(3.15) 


where their discreteness is implicit. 

The next section will discuss the inherent problems in using 
the DFT to approximate the Fourier integral transform and the complex 
Fourier series. 

3.2 The Discrete Fourier Transform 

When the DFT was introduced In the previous section, it was 
said to approximate the Fourier Integral transform and the complex 
Fourier series. The Interested reader not familiar with the derivation 
of the DFT and the nature of the approximations made is referred to 
several books which have been written on the subject (Ref. 10, 4). The 
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discussion of the DFT Included in this paper will he limited to the 
selection of sampling parameters which will avoid frequency and 
time-domain aliasing. 

V/hen utilizing the DFT, there are many sampling parameters 
which may be varied, including the sampling frequency, Nyquist frequen- 
cy, frequency resolution, length of the time record, time resolution, 
and the number of digital samples. However, using the relationships 
given in Appendix B, the specification of any two of these parameters 
will uniquely define the remaining parameters. Therefore, the two 
parameters will be chosen such that aliasing does not occur. 

In order to avoid frequency-domain aliasing, the Nyquist 
frequency Is selected in accordance with the following guidelines: 

* For band-limited functions - Select f N such that It Is 
greater than or equal to the highest frequency contained In 
the exponentially-windowed forcing function. 

* For functions whose spectrum is not band limited - Examine a 
typical system transfer function and determine the 
frequency, f , where the function effectively decays to 
zero, (i.e. H (f>f c ) = 0) . Since (X (f>f c )} = {0} 
in accordance with Eq. (3.15), a value for f^ is then 
selected such that the exponentially-windowed force spectrum 
below f is represented accurately. Example 1 in Chapter 
7 demonstrates this procedure. 
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Once the Nyquist frequency has been determined, the resolution 
in the time domain may then be found from 

T max = ^ (3 ‘ 16) 

The time resolution is denoted as a maximum resolution since it may be 
desired to select a smaller value of T In order that the system 
response may be obtained with a greater resolution. 

Now, In order to avoid time-domain aliasing, the length of the 
time record, T , will be chosen such that both the exponentially- 

A 

windowed excitation, {F(t) } , and exponentially-windowed response, 

A 

{x} , have 'effectively decayed to zero' before the end of the time 
record. The author has found that time-domain aliasing will not be 
significant If 'effectively zero' Is said to mean that the response at 
the end of the time record Is at least two orders of magnitude less 
than the response at the beginning of the time record. 

The remaining sampling parameters may now be determined by 
using the equations In Appendix B In conjunction with the values of f N 
and T q selected to prevent aliasing. The most important of these 
parameters is the number of digital samples, K , since It was seen in 
Eq. (3.2) that the value of K determines how many sets of frequency- 
domain equations must be solved. Thus, the minimum value of K may be 
determined from 


''n.tn ■ 2 Vn 


(3.17) 
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Once a suitable Integer value of K has been chosen, either T Q or 
f N must he recalculated since only two parameters may be independently 
specified. 

The DFT Is often implemented on the digital computer using the 
Fast Fourier Transform (FFT) algorithm (Ref. 10, 12) which has a higher 
computational efficiency than the DFT algorithm. Therefore, the 
following chapters will use DFT and FFT synonymously, since they differ 
only in implementation. 

It must be noted, however, that in a radix-2 implementation of 
the FFT, K is required to be an Integer power of two. Therefore, If 
the minimum value of K was determined to be 550, It would be neces- 
sary to use a K of 1024. For problems where n Is 'large', the cost 
of the solution of 474 extra sets of equations may be greater than the 
savings achieved by using the FFT algorithm. 

Before performing a DFT on a periodic signal, one additional 
topic requires discussion. That topic is leakage. In order to prevent 
leakage, there must be an integer number of fundamental periods of the 
signal within the time record. This condition can be difficult to 
satisfy if one is trying to keep f^ and T Q within a certain range 
to obtain a suitable value of K . A windowing function such as the 
Hanning or Tukey window may be applied to the function in an attempt to 
reduce leakage (Ref. 4, 10). 

In summary, the selection of the sampling parameters is very 
important when using the DFT. The following section will discuss the 
selection of the convergence factor . and its influence on the problem of 
time-domain aliasing. 
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3.3 The Selection of a Convergence Factor 

To prevent time-domain aliasing it Is required that both 

A A 

{f(t)} and {x} decay to 'effectively zero' by the end of the time 
record, or window as it is often called. In the previous section it 
was noted that T Q could be chosen such that this criterion would be 
satisfied for any value of a . However, it is often more efficient to 
chose T q to be the maximum time at which the response is desired. 
Since this value of T Q Is often less than the T q required to 
prevent time-domain aliasing for an arbitrary a , the value of K 
will be reduced in accordance with Eq. (3.17). 

Therefore, If T q Is selected without regard for the aliasing 
problem, the value of the convergence factor must be selected such that 

A A 

the functions {f(t)> and {x} decay to 'effectively zero' by the end 
of the window. An examination of Eqs. (2.26) and (2.31) reveals that 
an increased value of the convergence factor will cause the functions 
to decay faster and will further reduce the chances of time-domain 
aliasing. Therefore, It appears that the convergence factor should be 
chosen to be 'very large.' 

Although there is no limit on the maximum value of the conver- 
gence factor from an analytical standpoint, there is a limit Imposed 
from a computational standpoint. The problem occurs when large values 
of a are used such that the order of magnitude of the functions at 
opposite ends of the window are vastly different. This large variation 
In magnitude causes the computational model to be ill-conditioned in a 
finite word length computer, and therefore causes roundoff errors. 
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In a limited number of tests, this author concluded that the 
convergence function should not cause the original function to change 
more than two or three orders of magnitude within the window. A good 
rule of thumb for the maximum value of the convergence factor was found 
to be 


_ -ftn(10“ 2 /e 0 °) 


max 


(3.18) 


where a Q is the exponential order of the forcing function. This 
convergence factor causes the function to change two orders of magni- 
tude within the window. 

In summary, the value of the convergence factor must be greater 
than zero as determined in Chapter 2, and it must be large enough to 
prevent time-domain aliasing, while remaining small enough to prevent 
Ill-conditioning of the computational model. 

This concludes the discussion of parameter selection to obtain 
an efficient and well -conditioned frequency-domain model. The 
remainder of this chapter will give a synopsis of frequency-domain 
analysis and will discuss additional computational enhancements. 

3.4 A Synopsis of Frequency Domain Analysis 

At this point, the frequency-domain solution to steady-state 
and transient problems may be outlined as shown in Figure 3.1. The 
outline emphasizes the fact that [G^] , an n x n complex matrix, 
must be formed and inverted for each of the K discrete values of s . 
For a general loading condition where K=128 or 256. These 
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inversions may drive the cost of the analysis prohibitively high. This 
Is the most fundamental computational problem to overcome when using 
frequency-domain analysis techniques. 

In an attempt to bring the cost of the analysis down, the 
following suggestions to improve computational efficiency will be 
considered: 


* The reduction of K 

* The formation of [H^] directly (i.e. without calculating 
[G^] and Inverting it to obtain [H^]) 

* The movement of as many calculations as possible outside of 
the frequency loop 

* The reduction of the size of the matrices within the 
frequency loop 

The first suggestion may be evaluated by studying the prop- 
erties of the DFT of a real function. The study reveals that when a 
real function Is transformed using the DFT, the real part of the 
frequency spectrum has even symmetry about the Nyquist frequency, and 
the imaginary part of the spectrum has odd symmetry about the Nyquist 
frequency. Therefore, only the first K' points of the response 
function need to be calculated, where 



1 


K' 


( 3 . 19 ) 
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The remainder of the response spectrum may be determined from 

{X k ) = {X K _ k } * k = K'+l,. . . ,K-1 (3.20) 

* 

where {X K _ k ) is the complex conjugate of ( x K _ k > . Therefore, the 
effective K has been reduced to approximately half of its original 
value. 

The second and third suggestions will be evaluated in the 
following section, and the fourth suggestion will be evaluated in 
Section 3.6. 

3.5 Forming the Dynamic Flexibility Matrix Directly 

Forming the dynamic flexibility matrix directly can lead to a 
more efficient solution by allowing more computations- to be placed 
outside the frequency loop. It also eliminates the inversion of the 
complex dynamic stiffness matrix. While the procedure described in 
this section still requires a matrix Inversion, the equations forming 
the matrix are uncoupled yielding a diagonal matrix whose inversion is 
trivial. 

The direct formulation of the dynamic elasticity matrix will 
begin by writing the elgenproblem for the undamped structures, 

[-w 2 [m] + [k]] (x(w)} = {0} (3.21) 

A nontrivial solution of Eq. (3.21) may be obtained by setting 


det j-w 2 [m] + [k]| = 0 


(3.22) 
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and finding the roots. The roots obtained are the n eigenvalues of 
the system of equations in Eq. (3.21), which have n eigenvectors 
associated with them. The eigenvectors, <t> , are stored columnwise in 
a matrix which is denoted as $ . 

01 - (3.23) 

Equation (3.7) can now be transformed from physical space to 
modal space by using the transformation 

(\) a [*] (n k > (3.24) 

Substituting this into an expanded form of Eq. (3.7) and premultiplying 
by [*] T . the equations of motion in modal coordinates can be written 
as 

[CM] s£ + [C] s k + [K]] (n k > = l>? T {F k ) (3.25) 

where [M] and [Kl are n x n diagonal matrices commonly denoted as 
the generalized mass and stiffness matrices. The n x n matrix [C] 
will only be a diagonal matrix for special types of damping, where It 
can be shown that 

[c] or 1 [k] = [k] or 1 [c] (3.26) 

(Ref. 13). 

If Eq. (3.26) is satisfied, the equations of motion are un- 
coupled for each modal degree of freedom. At this point, Eq. (3.25) 
may be written in simplified form as 
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[A k ] {n k } = r*3 T (F k ) (3.27) 

where 

^ A k" = 1=l!n * M 1 s k + C 1 s k + K i ) ( 3 * 28 > 

Equation (3.27) may then be solved for {n k ) using the following 
equation 

(n k ) • CA^' 1 [*] T {F k > (3.29) 

However, since [A k l Is a diagonal matrix, its Inverse consists simply 
of the reciprocal of the diagonal terms. 

The equations may then be transformed back to physical space 
using the transformation In Eq. (3.24) 


where 

[D k ] 


<X k } = [*] [D k ] 0] T {F k J 
- [A k ] - 1 - ,2]** ( ! ) 


(3.30) 

(3.31) 


Upon comparing Eqs. (3.11) and (3.30), It Is realized that 
[H k ] may be calculated directly by 

CH k ] = [»] CD k ] of (3.32) 


Since the elgensolutlon and the formulation of the generalized 
mass, stiffness, and damping matrices are clearly frequency-indepen- 
dent; those operations may be placed outside of the frequency loop. It 
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should be noted, however, that the modal formulation of [H^] did not 
come without considerable expense. The major disadvantages of the 
modal formulation are: 

* An n x n elgensolution must be performed. 

* The class of problems has been restricted to those whose 
damping matrix is diagonalized by the eigenvectors. 

However, the advantage of the modal formulation is: 

* An n x n inverse is no longer required at each discrete 
frequency. 

These advantages and disadvantages must be weighed against each 
other for a particular problem. It is clear that If the equations need 
to be evaluated at only a single discrete frequency, the single inverse 
of [G^l would be more efficient than the eigensolution and subsequent 
matrix multiplications. 

3.6 Reducing the Problem 

The next suggestion to improve computational efficiency Is the 
reduction of the size of the matrices in the frequency loop. Recall 
that n is the number of degrees of freedom of the finite element 
model, and it cannot be reduced without reducing the accuracy of the 
solution. However, if the system response is desired at only a limited 
number of points, the number of system equations can be reduced without 
reducing the accuracy of the solution. 
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The reduction process Is initiated by partitioning the response 
vector into those degrees of freedom which are desired, R , and those 
which are not desired, I . The R-set will be denoted as the relevant 
set, and the I-set will be denoted as the irrelevant set (i.e. Irrele- 
vant in the sense that the response is not desired at those degrees of 
freedom). 

The partitioning yields 


{X} = 



(3.33) 


Equation (3.15) may now be partitioned and written as 


OC 

X 


h rr h ri 


V 


a 




X I 

> . 


h ir h ii 


F I 

» . 


(3.34) 


The response of the relevant degrees of freedom may be obtained by 
expanding the top equation, and writing it as 


{X R } = [H rr ] (F r > + tH RI ] (F r > 


(3.35) 


The problem may be simplified further by requiring all forces to be In 
the R-set. This will cause the {Fj} vector to be a null vector, and 
will allow Eq. (3.35) to be written as 


{Xr> = 

C H rr] {F r > 

(3.36) 

{X} = 

[H] {F} 

(3.37) 
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or 

[G] {X} = (F> (3.38) 

where the (") notation Indicates a matrix which has been reduced to 
its relevant degrees of freedom. The relevant part of the dynamic 
elasticity matrix, often called the dynamic elasticity transfer matrix, 
is a submatrix of the dynamic elasticity matrix, whereas the relevant 
part of the dynamic stiffness matrix, or dynamic stiffness transfer 
matrix, is not a submatrix of the dynamic stiffness matrix. 

For problems where the number of relevant degrees of freedom, 
p, is much less than the total number of system degrees of freedom, 
n , the use of Eq. (3.37) will be considerably more efficient than the 
use of Eq. (3.11) . The efficiency would be even higher, if [H] 
could be found directly (l.e. without forming [H] and removing the 
relevant portion). 

Geering presents the direct formulation of [H] which utilizes 
a projection scheme (Ref 3). However, by expanding Eq. (3.32), it may 
be shown that each element in the [H] , as well as [Hi , matrix may 
be obtained from 

H 1j = ^ *1 r *jr D r (3,39) 

where D r is the r^ diagonal element of [D^] defined in Eq. 
(3.31) and $.j r (or <j>j r ) represents the value at the i th (or j th ) 
row, and r th column of the modal matrix [♦] . 

It should be noted at this point that Eq. (3.37) yields an 
exact solution to Eq. (3.15) at the relevant degrees of freedom. It is 
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also an exact solution to £q. (2.1) if the forward transform of the 
dynamic load and the inverse transform of the dynamic response spectrum 
are exact. There Is, however, an additional reduction which can be 
made with an acceptable loss of accuracy. 

3.7 Further Reduction of the Problem 

The final reduction scheme presented in this chapter is the 
well-known modal approximation (Ref. 14, 15). It is important to note 
that this Is an approximation to the exact solution given by Eq. 
(3.37). The loss of accuracy, however, is generally considered to be 
acceptable. 

Modal approximation occurs when only N of the n system 
modes are retained In the modal expansion of the displacement field 
(where N < n). Equation (3.39) may then be written as 

N 

H.. = Z 4. *. D (3.40) 

ij ir y 1r r v ' 

r-i 

It Is Important to remember that if a system is reduced to p 
relevant degrees of freedom, then a minimum of p modes must be 
retained in the expansion to guarantee the invertlbllity of [H] . 

If a modal approximation is used, it is suggested that the 
residual modes also be included in the analysis In order to reduce the 
loss of accuracy (Ref. 14, 16). 

3.8 A Synopsis of Efficient Frequency-Domain Analysis 

At this point, the computational improvements may be incor- 
porated into the outline of Figure 3.1. Recall, however, that these 
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computational improvements are only improvements for a problem which 
must be evaluated at a large number of discrete frequencies, s^ . The 
new outline Is shown in Figure 3.2 

It is clear that the improvements will help reduce the cost of 
many analyses. The greatest computational concern now is the poten- 
tially large eigenvalue problem which must be solved In Step 2. The 
following chapter on substructure coupling will address this problem. 
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Determine the spectrum of the exponentially- 
windowed forcing functions at K discrete 
frequencies (Eq. (3.9)) 


2. Execute the following frequency loop K 
(k Index): 

a. Form [G^] (Eq. (3.8)). 

b. Invert [G^] to obtain [H k ] 

c. Solve for {X k > (Eq. (3.11)) 

- ★ 

3. Inverse transform (X) to obtain {x} 

★ 

the exponential window to obtain (x> 


times 


and remove 
(Eq. (3.13)) 


Figure 3.1 Outline for frequency-domain analysis 



1. Detemine the spectrum of the relevant exponentially- 
windowed forcing functions, {F> , at K discrete 
frequencies (Eq. (3.9)). 

2. Perform an eigenvalue analysis of the structure in 
order to obtain the modal matrix [*] . 

3. Form the [M] , [C] , and [K] matrices. 

4. Execute the following frequency loop K' times 
(k Index): 

a. Form [H^] (Eq. (3.38)) 

b. Solve for {X k > (Eq. (3.37)) 

5. Inverse transform {X} to obtain the relevant part 

* it 

of {x> and remove the exponential window to obtain 
the relevant part of {x}* (Eq. (3.13)). 


Figure 3.2 Outline of efficient frequency-domain analysis. 



Chapter 4 

SUBSTRUCTURE COUPLING IN THE FREQUENCY DOMAIN 

4.1 Equations of Motion for the Substructure 

Consider a superstructure 1 that has been separated Into two or 
more substructures. For Illustration purposes, consider the super- 
structure and Its substructures shown in Figure 4.1. 

Note that each substructure retains the same boundary condi- 
tions and external forces which it had as a part of the superstructure. 
In addition, each substructure has coupling forces at its Interfaces 
with other substructures. These coupling forces represent the forces 
transmitted by adjacent substructures. 

The equations of motion may be written Independently for each 
substructure as 


and 


[G 1 ] (X 1 ) = (F 1 > 
[G*] {X*} = {F*> 


(4.1) 

(4.2) 


The set of physical degrees of freedom for each substructure will now 
be divided into a set of interface or juncture degrees of freedom, J , 
and a set of external degrees of freedom, E , which will encompass 
the remaining active degrees of freedom of the substructure. The 
substructure equations of motion may be written in partitioned form as 


1 This paper will denote the entire structural system either as the 
superstructure, or as substructure 0 . 
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Substructure 0 (Superstructure) 



Substructure 1 Substructure 2 


Figure 4.1 Typical substructure model . 
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(4.3) 


(4.4) 


where the following notation applies to the i* 11 substructure 


Xj Interface DOFs 

External DOFs 

f£ Interface forces due to adjacent substructures 

Fj External forces at Interface DOFs 

F^ External forces at external DOFs 


The frequency-domain equations for each substructure may now be 
coupled together to form the superstructure eouations of motion. This 
coupling Is accomplished by enforcing the displacement and force 
compatibility at the substructure Interfaces. 

The displacement compatibility at the interface may be ex- 
pressed as 


{x j } s {x j } = {X °> (4.5) 

Equation (4.5) ensures that the Interface displacements of adjacent 
substructures are identical and are equivalent to the displacement of 
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the same degrees of freedom on substructure Q (i.e. the superstruc- 
ture). Force compatibility at the Interface may he written as 


< F C> ■ -< F c } 


(4.6) 


At this point, Eqs. (4.3) and (4.4) may be reordered and 
combined to yield 


g ee 

g ej 




rx 1 

A E 


\4 1 

•Se 

g L 




X 1 

A J 


F i tF i 



r 2 

G 0J 

r 2 

g je 


X 2 

A J 


f j + f I 



r 2 

g ej 

r 2 

g ee 


X 2 

A E 


ii 1 


By adding the second and third equations together, and enforcing the 
displacement compatibility condition, the system may be reduced to 
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(4.8) 


Now the force compatibility condition may be enforced to yield the 
equations in their final form (Ref. 3). 
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G 
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JE 



(G JJ + G Jj) 






(4.9) 


This process of enforcing compatibility of the forces and 
displacements at the substructure Interfaces is analogous to the direct 
stiffness approach used In static analysis. The only difference In 
this dynamic formulation Is that the stiffness matrix, displacement 
vector, and force vector are not constants, but are, instead, functions 
of frequency. 

Extending the analogy of the direct stiffness method, the 
element stiffness matrix, or the substructure dynamic stiffness matrix, 
will be defined as 

[G 1 ] = [m 1 ] sj; + fc 1 ] s k + He 1 ] (4.10) 

where [m 1 ] , Tc 1 ] , and [k ] are the free-interface mass, damping 
and stiffness matrices of the 1 substructure . Similarly, the 

element elasticity matrix, or the substructure dynamic elasticity 
matrix, will be defined as 

•x 

Recall from Chapter 2 that only active degrees of freedom are 
included In the equations of this paper. 
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CHji - r*'i [dJ] iVi t (4.1D 

where [4 1 ] Is the substructure free-lnterface modal matrix, and (0.1 
may be found by substituting the substructure's generalized mass, 
damping, and stiffness matrices Into Eq. (3.31). As demonstrated in 
Chapter 3, Eq. (4.11) may be expanded and each term of the dynamic 
elasticity matrix may be expressed as 

’ j 1 *l*3r D r < 4 - 12 > 

The number of equations represented by Eq. (4.9) may be cal- 
culated by adding together the degrees of freedom of each substructure, 
m 1 , and subtracting the total number of Interface degrees of freedom, 

Pl ° • or 

M i 

n = i m - p; (4.13) 

1=1 1 

where M Is the total number of substructures. 

At this point. It should be noted that Eq. (4.9) is simply a 
partitioned representation of Eq. (3.14). Therefore, the same solution 
procedures and modifications described In Chapter 3 are applicable here 
as well. 


4.2 A Synopsis of Substructure Coupling 

As in Chapter 3, the above procedure has been outlined In order 
to determine where the method could be improved upon. An evaluation of 
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the outline shown in Figure 4.? again reveals that a reduction in the 
amount of computation within the frequency loop will yield a more 
efficient solution to the problem. 

The same computational improvements made at the system level 
apply equally well here, at the substructure level. However, the 
interface degrees of freedom must now be Included in the relevant set 
In order to provide the coupling to adjacent substructures. Therefore, 
the relevant set will now consist of Interface degrees of freedom as 
well as forced degrees of freedom and other degrees of freedom where 
the response is desired. 

The substructure coupling outline with the improvements from 
Chapter 3 Is shown in Figure 4.3. An examination of the outline shows 
that the use of substructures has required additional matrices to be 
inverted in step 3. a. 2. Since the number of calculations required to 
perform a matrix inversion is on the order of b 3 (where b Is the 
order of the matrix), it is essential that the size of those matrices 
be reduced. 

• * • • 

Recall that [Hj Is a p x p matrix, where p 1 Is the 

number of relevant degrees of freedom of the 1 th substructure. It has 
already been stated that interface degrees of freedom must be Included 
in the relevant set. Therefore, in order to minimize the number of 
degrees of freedom In the R-set, one must transform all loads to the 
interface and request the response at only the interface degrees of 
freedom. 

However, requesting the response at the interface degrees of 
freedom only, is not an acceptable solution, unless the response can be 


1. Peteraine the spectrum of the exponentially- 
windowed forcing functions at K discrete 
frequencies (Eq. (3.9)). 

2. Execute the following frequency loop K times 
(k index): 

a. Execute the following substructure loop 
M times (i index): 

1) Form [Gj] Eq. (4.10) 

2) Using direct stiffness methods, sum the 

matrices to form [g£] , the 
superstructure dynamic stiffness matrix. 

b. Invert [g£] to obtain [h£] 

c. Solve for {X k > (Eq. (3.11)) 

* * 

3. Inverse transform (X> to obtain {x} and 

★ 

remove the exponential window to obtain {x> 

(Eq. (3.13)). 


Figure 4.2 


Outline of substructure coupling in the 
frequency domain. 


1. Determine the spectrum of the relevant 
exponentially-windowed forcing functions at 
K discrete frequencies (Eq. (3.9)). 

2. Execute the following substructure loop M 
times (i index): 

a. Perform a free-interface eigenvalue 
analysis of each substructure in order 
to obtain fV] . 

b. Form the [M^3 , [C^l , and [K 1 ] 
matrices. 

3. Execute the following frequency loop K' 
times (k Index): 

a. Execute the following substructure loop 
M times (i Index): 

1) Form [hJ] (Eq. (4.12)) 

2) Invert [h£] to obtain [G^] 

3) Using direct methods, sum the TG^] 
matrices to form [g£] , the super- 
structure dynamic stiffness matrix. 

b. Invert [g£] to obtain [h£] 

c. Solve for (X k ) (Eq. (3.11)) 

A if 

4. Inverse transform (X) to obtain lx) and 

★ 

remove the exponential window to obtain (x) 
(Eq. (3.13)) 


Figure 4.3 Outline of efficient substructure coupling. 
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recovered subsequently at the irrelevant degrees of freedom. This 
topic is discussed in the following section. Reducing the loads to 
Interface degrees of freedom is then discussed In section 4.4. 

4.3 Recovery of the Response at Irrelevant Degrees of Freedom 

Recovery of the irrelevant degrees of freedom may be accom- 
plished by first writing the substructure eouations of motion parti- 
tioned Into their relevant and Irrelevant sets 

(4.14) 

The bottom equation may then be expanded and solved for the response at 
the irrelevant degrees of freedom. 

[Gjr] (X r > + [Gjj] {Xj} = (Fj> (4.15) 

{Xj} - [Gjj ]' 1 ({Fj> - [Gj R ] {X r }) (4.16) 

For the special case where all loaded degrees of freedom have 
been placed in the R-set, (Fj > Is a null vector and the response of 
the irrelevant degrees of freedom may be determined using the simpli- 
fied equation 

(Xj) = -[[Gji]’ 1 [G ir ]] (X r > (4.17) 

where the {X R } vector was determined previously by the evaluation of 
the superstructure response. 
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For many problems, this method of calculating response is more 
efficient than Including the desired response points in the system 
relevant set. As an example, consider two 50 degrees-of -freedom 
substructures attached at 10 Interface degrees of freedom. If the 
response is desired at each of the 90 system degrees of freedom, and 
all degrees of freedom were retained in the relevant set, the inversion 
required In step 3.b of Figure 4.3 would require on the order of 90 3 or 
729,000 operations. However, if only the Interface degrees of freedom 
were retained as relevant degrees of freedom, the inversion in step 3.b 
would be on the order of 10 3 or 1,000 operations. Of course, the 
subsequent postprocessing required in order to obtain the response at 
the remaining 80 degrees of freedom would require the Inversion of two 
40 degree of freedom matrices increasing the cost by another 128,000 
operations. However, the total processing cost of the latter method Is 
still considerably less than including all degrees of freedom in the 
relevant set. However, if the cost of determining the response were to 
be examined for the same substructures attached at 40 degrees of 
freedom, one would find that is was more efficient to retain all 
degrees of freedom in the relevant set. 

Although the efficiency of the two methods is problem depen- 
dent, the more efficient one may be determined a priori by examining 
the cost of the inversions required by each of the two methods. 

4.4 Transfer of Substructure Loads to the Interface Degrees of 

Freedom 

This final reduction scheme in Chapter 4 will allow all of the 
substructure's external loads to be transformed into equivalent inter- 
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face loads. This reduction is very valuable for the class of problems 
with a large number of external forces (e.p. acoustics, distributed 
loads, etc.). 

This procedure will require another matrix inversion, therefore 
the system should be reduced so that only interface and loaded degrees 
of freedom are in the relevant set prior to beginning this reduction. 
The interface degrees of freedom will again *orm the J-set, and the 
remaining degrees of freedom will form the E-set. The substructure 
equation of motion may then be partitioned and written as 

(4.18) 

Since the interface degrees of freedom will make up the R-set, 
and the external degrees of freedom will make up the I-set, Eq. (4.18) 
may be written as 

(4.19) 

The bottom equation may be expanded and solved for (X j } . 

fXj} « [Gjj]” 1 ({Fj> - [G ir ] (X r » (4.20) 

At this point, the top equation is expanded using the substitution for 
{Xj > derived above. 
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CGrrD (X r ) + [G ri 1 [Sj,]* 1 «Fjt - fG !R ] {X R » ■= (F r ) (4.21) 

Simplifying the notation, 

[G] (X) = (F) (4.22) 

where 

{X} = (X R ) (4.23) 

and 

[G] s [G RR ] ' [G R j] ^ G II^ * [Gj R ] (4*24) 

and 

(F) = (F r > - [G ri ] [Gjj]" 1 {Fj> (4.25) 

When this reduction scheme is applied to static (l.e. constant) 
matrices. It Is called 'static condensation.' Since the reduction 
shown above uses dynamic (l.e. not constant) matrices. It is appropri- 
ate to call it a 'dynamic condensation.' 

It Is Interesting to note the similarity between Eq. (4.22) and 
Eq. (3.38). The fact is, they are not only similar, but they are 
Identical. Therefore, the direct formulation of [H] may also be 
viewed as a dynamic condensation of [G] to obtain [6] , and an 
Inversion to obtain [H] . This fact will prove useful when deter- 
mining the natural frequencies of a system which has been reduced to 
its relevant degrees of freedom. 



Chapter 5 

DETERMINING THE NATURAL FREQUENCIES OF UNDAMPED STRUCTURES 

5.1 The Undamped Free Vibration Problem 

The natural frequencies of an undamped structure may be found 
by examining the free-vibration equations of motion for that structure, 

[m] (x) + [kl {x} = {0} (5.1) 

This equation may be transformed to the frequency domain using the 
methods described In Chapter 2. The free-vibration problem may then be 
described by the continuous frequency-domain equation 

[G(s)] (X(s ) } = {0} (5.2) 

where 

[G(s)1 = [m] s 2 + [k] (5.3) 

It has been shown that a nontrivial solution to Eq. (5.2) must satisfy 
the relationship 

det |[G(s)]| = 0 (5.4) 

where |*| is the determinant operator. 

An expansion of Eq. (5.4) will yield an n th -order polynomial 
equation in the complex frequency variable, s 2 . It has been shown 
(Ref. 17) that if [kl and [m] are real, symmetric matrices, where [ml 
is positive definite and [k] is either positive definite or 
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positive semi -definite, then the roots of Eq. (5.4) will be purely 
Imaginary. Thus the substitution 

s = joj (5.5) 

may be used In Eqs. (5.2) and (5.3). However, since s (and therefore 
jo>) appears In Eq. (5.3) as a squared term only, the equations become 
purely real and may be written In terms of the real variable u . Thus 

£6(«)] (XU)} - {0} (5.6) 

where 

[G(w)] = -a, 2 [m] + [k] (5.7) 

Finally, Eq. (5.4) may be written as 

det |[G(<a)]| ■ 0 (5.8) 

Either Eq. (5.6) or Eq. (5.8) may be used as the basis for obtaining 
the natural frequencies of the system as will be shown in the following 
discussion. 

The first approach Is to expand Eq. (5.6) and write It In the 
form of the generalized eigenvalue problem, 

Ck] M = [m] l>] (5.9) 

The square roots of the n real eigenvalues of this equation are the 
natural frequencies of the system. There are many algorithms which may 
be used to solve this generalized eigenvalue problem (Ref. 17, 18). 
Although a discussion of these algorithms is beyond the scope 
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of this thesis, the problems associated with obtaining the [m] and 
[k] matrices will be discussed. 

The second approach will address the solution of Eq. (5.8). 
The algorithms discussed in this paper do not employ the [m] and [k] 
matrices directly. Instead, the algorithms evaluate the determinant of 
[G(oj)] at discrete numerical values of u . The n positive real 
values of <o which satisfy Eq. (5.8) are the natural frequencies of 
the system. 

The next section will discuss how the [m] and [k] matrices 
may be recovered for unreduced frequency-domain problems, and Section 
5.3 will explain why the recovery is not practical for a reduced 
problem. The remaining sections will then discuss the second approach 
mentioned above. 

5.2 Recovering the System Mass and Stiffness Matrices 

When using the substructure coupling technique outlined in 
Section 4.1, the [m] and [k] matrices were not calculated for the 
superstructure. Only [Gl and [H] were known at the system level. 
It is reasonably simple, however, to recover the system [m] and [k] 
matrices for this class of frequency-domain problems where the system 
has not been reduced. 

The first step taken to recover the [m] and [k] matrices is 
to write the equations which define [G] for the undamped system at 
two distinct frequencies u>j and Ug , 

[G(uj)1 = -uj 2 [m] + [k] 


(5.10) 
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CGCto^) 1 = -to 2 2 Tm] + [k] (5.11) 

Since the mass and stiffness matrices are invariant, they are identical 
in each of the equations above. Therefore, Eq. (5.11) may be solved 
for [k] , 


[k] = [Gfa^)] + oj^ 2 [m] (5.12) 

which may then be substituted into Eq. (5.10) to yield 

[G(a)j)] = -ajj 2 [m] + [G(a> 2 )] + ^ ^ (5.13) 

Finally, Eq. (5.13) may be solved for [m] and written as 

[ml » [6(^)1 - [G(o» 2 )1 (5.14) 

(<d 2 2 - 0)j 2 ) 

The [m] and [k] matrices obtained from Eqs. (5.14) and 
(5.12) may then be used In any suitable eigensolution to obtain the 
system natural frequencies. It must be noted that although Eqs. (5.12) 
and (5.14) are algorithmically correct, their stability is not guaran- 
teed for any arbitrary pair of w which might be chosen. 

5.3 Effects of Reduction on the Natural Frequencies 

Up to this point, methods of obtaining the natural frequencies 
have only been considered for systems where all degrees of freedom have 
been In the relevant set. However, to utilize frequency-domain sub- 
structuring techniques effectively, only a minimum nuntoer of degrees of 
freedom are retained in the relevant set. 
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In Chapter 3, it was determined that the system could be 
reduced by keeping only the rows and columns of the dynamic elasticity 
matrix which correspond to the relevant degrees of freedom. It was 
also shown in Chapter 3 that this reduction scheme is equivalent to a 
dynamic condensation of the dynamic stiffness matrix. 

Since each of the reduction techniques is exact (i.e. no 
approximations are made), it is reasonable to assume that all of the 
dynamic characteristics are retained in the dynamic stiffness transfer 
matrix. 

Consider the expansion of [G] , given by Eq. (4.24), in terms 
of the relevant and Irrelevant partitions of the [m] and [k] 
matrices, 

[G(o, k )1 = (-<«>k 2 ^RR^ * I^rr]) ~ (-w^ 2 f- m Rj- + f k RjJ) * 

(5.15) 

(-u > k 2 C m ii] + C k ii ^ _1 (‘“k 2 ^ m IR^ + ^ k IR^ 

Although each of the partitions which make up the system [m] 
and [k] matrices are included in the equations, it is not immediately 
clear how they could be recovered for use in Eq. (5.9). Furthermore, 
even if the system [m] and [k] matrices were recovered, an n x n 
elgensolution would then be required to obtain the system natural 
frequencies and the computational advantages of substructuring would be 
lost. 

Therefore, the remainder of the chapter will discuss an alter- 
nate solution procedure for determining system natural frequencies 
which can be extended to handle reduced systems. 
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5.4 Determining the System Natural Frequencies 

Consider Eq. (5.8) upon which this procedure for determining 
natural frequencies will be based. The expansion of this equation 
yields an n^-order polynomial in u 2 . The n positive roots of this 
polynomial (which is often called the characteristic equation) 
correspond to the n natural frequencies of the system. A plot of a 
typical characteristic equation is shown in Figure 5.1. 

One of the most popular methods for determining the natural 
frequencies of a system described by Its characteristic equation is the 
determinant search method (Ref. 18). A pure determinant search method 
searches for zero crossings of the characteristic polynomial by noting 
where the determinant of [G] changes sign. Once a root Is bracketed 
(l.e. the determinant is positive on one side of the root and negative 
on the other side). It is relatively easy to determine its value using 
a bisection scheme. 

There are, however, problems which arise when using a pure 
determinant search method. They are: 

It is Impossible to locate repeated natural frequencies. 

It is difficult to locate closely-spaced natural 
frequencies. 

It is not possible to determine how many natural 
frequencies lie within a given frequency range. 

From a computational point of view, the method is not very 
efficient either. The outline for the procedure, shown in Figure 5.2, 
reveals that the incremental search procedure requires a large number 
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Figure 5.2 Outline of the pure determinant search method. 
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of determinants to be evaluated for problems where the search bandwidth 
Is large, but the frequency increment Is small in order to be able to 
bracket closely-spaced natural frequencies. Therefore, the key to 
reducing the number of determinant evaluations is to come up with a 
better technique for bracketing the roots. 

The following section discusses an enhancement to the determi- 
nant search method which solves both the bracketing problem, and the 
problem the pure method has with locating all of the roots of the 
characteristic equation. 

•• 

5.5 Utilizing Sturm Sequence Checks 

One of the most popular means of enhancing the determinant 

•• 

search method Is by incorporating Sturm sequence checks into the 

•« 

determinant search algorithm. The Sturm sequence check allows each of 
the problems noted In the previous section to be overcome by deter- 
mining the number of natural frequencies below a specified frequency. 

•• 

This number, denoted as the Sturm sequence count, S , may then be 
used In conjunction with a bisection scheme to bracket each zero of the 
characteristic equation. 

•* 

A detailed discussion of the Sturm sequence check is beyond the 
scope of this thesis, and thus the Interested reader is referred to the 
work of Bathe and Wilson (Ref. 18). At this point. It. will suffice to 
know that the Sturm sequence count is equal to the number of negative 
elements on the diagonal of [D] , where [D] is obtained from the 
[L][D][L] T decomposition of [G] . (Note that the explicit frequency 
dependency notation has been dropped to simplify the notation.) Thus, 
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where 


n 

S = z c. 
1=1 1 


■ 1 

if 

d n 

< 0 

0 

if 

d ii 

2 0 


(5.16) 


where d^. represents the term on the diagonal of [D] . 

•• 

The determinant search method utilizing Sturm sequence checks 
outlined in Figure 5.3 will now be able to located all roots of the 
characteristic eouatlon, and It will do so with a much higher effi- 
ciency than before. 


5.6 Determining the Natural Frequencies of Reduced Systems 

Up to this point, the determinant search method has only been 
discussed for unreduced frequency-domain models. However, in order to 
take full advantage of the computational efficiency of using reduced 
models, the following extension of the determinant search method Is 
presented. 

First, consider the free vibration equations of motion for the 
undamped system reduced to Its relevant degrees of freedom 

[GT {X> = {0} (5.17) 


In general, a nontrivial solution to this problem must satisfy 


det | G ( = 0 


(5.18) 
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Figure 5.3 Outline for determining the natural frequencies 
utilizing Sturm sequence checks 
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An expansion of Eq. (5.17) and an examination of Figure 5.4 will reveal 
that the characteristic equation now has poles as well as zeros and is 
not a simple n th -order polynomial in a> 2 as before. This may be 
explained by examining Eq. (4.24), which reveals that the poles occur 
at frequencies where [Gjj] Is singular, and therefore its inverse 
approaches infinity causing the determinant of Tfi] to approach 
Infinity also. These frequencies may be thought of as the natural 
frequencies of the system which has all of its relevant degrees of 
freedom constrained. 

It Is obvious that a pure determinant search method would find 
this environment difficult to work in, since the characteristic equa- 
tion may now change sign via Infinity, as well as via zero. Further- 

•• 

more, the standard Sturm sequence check will not locate all of the 
roots either, due to reasons which are detailed In the following 
discussion. 

Since It Is the natural frequencies of structures composed of 
substructures which Is of primary interest in this thesis, consider a 
superstructure which has been reduced such that only the Interface 
degrees of freedom are In the relevant set. The equation of motion 
describing this system has been stated previously In Eq. (5.17). 

As In the previous section, the Pj° x Pj° dynamic stiffness 
transfer matrix, [G] , may be factored using an [L][D][L] T 
decomposition. By taking into account that the forcing function Is a 
null vector, the decomposed equation of motion may be simplified and 


written as 
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[U1 (X} = {0> (5.19) 

where 

[0] = [6] [[] T (5.20) 

and [0] Is an upper triangular matrix. 

•• 

The Stum seauence count of the reduced system may now be 
obtained by counting the number of negative elements on the diagonal of 
[5] . Thus 


where 



1 

If 

d ii 

< 0 

0 

if 

d 'n 

o 

Al 


(5.21) 


and d^. represents the i th diagonal element of [6] . It should be 

•• 

clear at this point, that the maximum attainable Stum sequence count 
for this reduced system is Pj° since there are not more than Pj° 
diagonal elements on the diagonal of a Pj° x Pj° matrix. 

•• 

In order to determine the significance of the Stum sequence 
count calculated in Eq. (5.21), write the equations of motion of the 
unreduced system partitioned into its relevant and irrelevant parts. 



(5.22) 



Equation (5.22) may now be decomposed and written in Its upper diagonal 
form, 
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(5.23) 


Since [U] is an upper triangular matrix, [U RI ] = [0] and the 

bottom equation may be written as 


[U RR ] {X R > = {0} (5.24) 

Recalling that {X R } = {X} , the comparison of Eqs. (5.19) and 

(5.24) reveals that 


[0] - [U RR ] (5.25) 

Therefore, the Sturm sequence count of the relevant degrees of freedom 
in Eq. (5.24) is equivalent to the Sturm sequence count of the reduced 
system in Eq. (5.19) (Ref. 6). Now, only the significance of the Sturm 
sequence count of the irrelevant degrees of freedom needs to be 
discussed. 

Consider the definition of the Sturm sequence count for the 
(n - Pj°) degrees of freedom in the Irrelevant set. 

The Sturm sequence count for the (n - Pj ) degrees 
of freedom Is equivalent to the total number of 
fixed-interface substructural natural frequencies 
which are less than the specified frequency. 
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Thus, the count for the irrelevant degrees of freedom may be written as 

cj (5.26) 

i 

if to", < a) 
r 

if J 2 a) 
r 

and represents the r 1 '" fixed- interface natural frequency of the 

f K A 

j substructure, and u is the freouency at which the count is to be 
determined. 

•• 

In summary, the effective Sturm sequence count, S , may be 
determined by 

M i 

S » S + r Si B (5.27) 

j=l 

This effective Sturm sequence count Is equivalent to the Sturm sequence 
count obtained by Eq. (5.16) for the unreduced system, and It may be 
used In conjunction with a bisection scheme to determine the natural 
frequencies of the superstructure. 

5.7 Conclusions 

This method of determining the natural frequencies of a super- 
structure using frequency-domain substructure models has both advan- 
tages and disadvantages over the time-domain methods. 


where 


S^ 

^SB 


fir^-p^ 

z 

r=l 


10 





70 


The advantages are: 

There are no large superstructure elgensolutions to 
perform. 

The method allows for the calculation of only a specified 
number of natural frequencies (e.g. the first 10 natural 
frequencies) or the natural frequencies within a given 
frequency range. 

and the primary disadvantage is: 

The fixed-interface natural frequencies must be determined 
for each of the substructures using an elgensolution or 
other suitable method. 

It must be remembered however, that the number of calculations 
involved In an elgensolution Is greater than n 3 . Therefore, the cost 
of several small elgensolutions plus the cost of the frequency-domain 
search described in this section is often less than that required to 
perform the superstructure eigensolution. As Illustrated in Section 
4.3, the cost of the analysis is dependent upon the number of Interface 
degrees of freedom and other problem-dependent parameters. 



Chapter 6 

DETERMINING THE NATURAL FREOUENCIES OF DAMPED STRUCTURES 

6.1 The Damped Free Vibration Problem 

The natural frequencies of damped structures may be found using 
a frequency-domain technique similar to the one discussed in Chapter 5. 
The development of the technique will follow a path similar to that in 
Chapter 5 and begins by first writing the free-vibration equations of 
motion for a damped structure In the time domain. 


[m] (x> + [c] (x) + [k] (x> = {0} (6.1) 

and in the frequency domain, 

[G(s)1 (X(s)> = {0} (6.2) 

where 

[G(s)J = [m] s 2 + [c] s + [k] (6.3) 

As in Chapter 5, a nontrivial solution to this equation must satisfy 

det |G(s) | - 0 (6.4) 

This equation may also be expanded to yield an n^-order 
polynomial equation In s 2 having 2n roots. However, with the 
inclusion of the damping matrix, the roots nay be complex as well as 
real . 

Either Eq. (6.2) or Eq. (6.4) may be used to obtain the damped 
natural frequencies of the system. However, the solution of Eq. (6.2) 
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is subject to the same problems of obtaining the system mass and 
stiffness matrices of reduced systems as seen in Chapter 5. Therefore, 
this chapter will concentrate on the solution of Eq. (6.4) which is 
more applicable to frequency-domain models. 

6.2 Determining the Damped Natural Frequencies of Reduced Systems 
Since it is the determination of natural frequencies of reduced 
systems which is of primary interest in this thesis, this discussion 
will begin by considering the reduced frequency-domain equations of 
motion for a damped system 

[G] {X} = {0} (6.5) 

Although this equation appears identical to Eq. (5.17), the Pj° x pj° 
dynamic stiffness transfer matrix is now complex Instead of real as it 
was in Chapter 5. Since the transfer matrix is complex, a general 
nontrivial solution must now satisfy 

det |G(s) | = 0 (6.6) 

where both the real and imaginary parts of the determinant are zero. 

This chapter will also employ a determinant search technique to 
solve Eq. (6.6). However, the search for the roots must now be made In 
the complex s-plane, whereas In Chapter 5, It was made only along the 
imaginary axis of the s-plane. 

The search for the 2n complex roots may be restricted to the 
positive Imaginary half-plane by noting that If a complex root is a 
solution to Eq. (6.6), then the complex conjugate of the root is also a 
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solution. Therefore, the search may be restricted to quadrants I and 
II of Figure 6.1. The search may be restricted further by noting that 
for conservative systems, the real part of the natural frequency must 
be less than zero, and therefore, the search may be restricted to 
quadrant II. 

A pure determinant search method for locating the roots of Eq. 
(6.4) in the complex plane is subject to the same problems with identi- 
fying repeated roots and closely-spaced roots as before. It Is also 
possible that the determinant will change signs via Infinity as before, 

as shown In Figure 6.2. Therefore an Index-checking scheme similar to 
•« 

the Sturm sequence checks used In Chapter 5 Is also desirable when 
searching for roots In the complex plane. 

•« 

6.3 Extending Sturm Sequence Checks to the Complex Plane 

*• 

An examination of the derivation of the Sturm sequence check 
reveals that one of the most fundamental properties upon which the 
method Is based Is the eigenvalue separation property of undamped 
systems. The property may be stated as (Ref. 18) 
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(6.7) 


In essence, the property says that applying an additional constraint to 
the system will cause the eigenvalues of the constrained system to fall 
between those of the system without the additional constraint. 
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Although the separation property was derived for an undamped 
system, it will be assumed that it can be extended to damped systems by 
the following 

I Re (sj k> )l s | Re {sj k+1> )| s |Re (s| k> )| ••• (6.8) 

and 

|Im(s^ k> )| s |Im <sj k+1> )| s |Im (s| k> )| ••• (6.9) 

where s r is the r th damped natural frequency (i.e. s r = a r - jay) . 

A mathematical proof of Eqs. (6.8) and (6.9) is beyond the scope of 
this thesis; therefore, the equations are presented without proof. 
However, an example will be presented in Chapter 7 which will demon- 
strate the applicability of the equations to one particular system. 

Recall that for the undamped system, the dynamic stiffness 
matrix [6] is real, and that the Sturm sequence check was used to 
locate the Imaginary parts of the roots (the real parts were known to 
be zero). Although [G] is complex for the damped system, Eq. (6.9) 
will permit the Sturm sequence check to be used to determine the 

Imaginary part of the root from the real part of the triangularized 

•• 

[G] matrix. Likewise, Eq. (6.8) will permit the Sturm sequence count 
of the Imaginary part of the triangularized [G] to be used to obtain 
the real part of the root. The above statements may be expressed as 


( 6 . 10 ) 
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where 


and 


where 


C 1 


1 1* Re (d..) < 0 
0 If Re (d..) 2 0 



f6.ll) 


1 if In (d^) < 0 
0 if Im (d i - ) * 0 


where «n is the 1 th term of the diagonal matrix resulting from the 

ClJ[D][L] T decomposition of [G] . 

Reduction of the system and substructure ng presents the same 

problems for the damped system as for the undamped. An argument 

similar to the one used In Section 5.6 may be used to derive the Sturm 

sequence counts for the damped system with substructures. For the sake 

of brevity, only the final equations are presented here. 

•• 

The Sturm sequence count which will be used to determine the 
imaginary part of the root may be written as 


I -I M T j 

S 1 = S 1 + z (sL) 

j=l 58 

where 



( 6 . 12 ) 


(6.13) 
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and 


1 if Re (d^.) < 0 
0 if Re (d if ) 2 0 



1 if Im (s r ) < 0 

0 if Im (s„) & 0 
r 


Similarly, the Sturm sequence count which will be used to 
real part of the root may be written 


where 


and 


R - D M P j 

S K = S K + Z (S" B ) 


j=l 


-R Pt - 

S K = Z l b 
i=l 


1 


1 if Im (d jf ) < 0 


0 if Im (d^) 2 0 


, nH-p 


3 


(S SB> J ' »r 


(6.14) 


determine the 


(6.15) 


(6.16) 


(6.17) 
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1 if Re (s r ) < 0 

0 if Re (s r ) * 0 

The next section describes an algorithm which has been imple- 
mented to demonstrate the applicability of Eqs. (6.8) and (6.9) in 
finding complex natural frequencies. The algorithm utilizes the Sturm 
sequence counts obtained in Eqs. (6.12) and (6.15). 

6.4 Bisection in the Complex Plane 

The damped natural frequencies may now be determined using the 

•• 

Sturm sequence indices given by Eqs. (6.12) and (6.15). However, a 

simple bisection scheme will no longer suffice for determining the 

roots, since the bisection of a plane is not unique. One relatively 

simple method of determining the roots is to combine a bilinear search 

algorithm with a bisection algorithm. 

The objective of the search is to find the value of s in the 

R I 

complex plane where both S and S change from (r-1) to (r) 

when looking for the r th root. This frequency, s , corresponds to 

the r tfl natural frequency of the system. Figure (6.3) has been 

•• 

Included to demonstrate how the Sturm sequence counts vary In the 
search quadrant of the complex plane. 

In short, as the bilinear search algorithm is looking for the 
r th root, it holds the real part of s constant and then uses bisec- 
tion to obtain an estimate of the imaginary real part of the r*^ 1 root. 
It then holds the imaginary part of s constant and uses bisection to 
obtain an estimate of the real part of the root. This procedure is 
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then repeated until both searches satisfy the tolerances specified. 
The value of s at which both tolerances are satisfied Is the r th 
damped natural frequency of the superstructure. 

6.5 Conclusions 

From a computational standpoint, the advantages and disadvan- 
tages of using frequency-domain techniques Instead of time-domain 
techniques to determine the damped natural frequencies of a system are 
similar to those listed in Section 5.7 for an undamped system. The 
only significant difference is that the damped natural frequencies are 
now required for each fixed interface substructure. 

From a practical standpoint, the primary concern with using 
this frequency-domain method is that Eqs. (6.8) and (6.9) have not been 
proven explicitly. The numerous examples which have been examined by 
the author have been successful. However, this Is not a proof that the 
equations are valid for all systems. 



Chapter 7 
EXAMPLES 


7.1 Introduction 

This chapter has been included to demonstrate some oT the 
frequency-domain techniques discussed In previous chapters. The 
examples do not detail the computations involved in the analysis of a 
problem, but instead they outline the major steps in each of the 
frequency domain analysis techniques and demonstrate the selection of 
proper analysis parameters. 

The following examples are Included: 

Example 1 - This example demonstrates how the 
transient response may be determined using frequency- 
domain techniques, and the importance of selecting 
the proper sampling parameters. 

Example 2 - This example demonstrates the coupling 
of substructures In the frequency domain and provides 
additional details on the selection of sampling 
parameters. 

Example 3 - This example demonstrates how the natural 


frequencies of undamped superstructures may be ob- 
tained using frequency domain substructure models. 
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Example 4 - This example demonstrates the applica- 
bility of Eqs. (6.8) and (6.9) for determining the 
damped natural frequencies of a superstructure using 
frequency domain substructure models. 


Example 1 

The first example will demonstrate the importance of careful 
selection of the sampling parameters as discussed in Chapter 3. The 
structure to be used In this example is the simple 1-DOF mass-spring 
oscillator 



with the following system parameters: 

n = 1.0 lbs 2 /1n * 31.42 rads/sec 

k = 987.0 Ib/in f x = 5.0 Hz 

The transient response of the system to two different dynamic loadings 
will be determined using the procedures outlined in Chapters 2 and 3. 


A. Selecting the Nyqulst Frequencies 

In order to demonstrate the selection of a proper Nyqulst 
frequency, the structure will be subjected to the cosine pulse loading. 


p(t) = 


0 t < 0 

987 cos(lOirt) 0 s t s 0.05 


0 


t > 0.05 
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A pulse such as this is not band-1 ini ted, and therefore the Nyquist 
frequency cannot be selected as the highest frequency in the 
exponentially-windowed pulse. Therefore, the alternate method of 
determining f N described in Chapter 3 will be used. 

Before proceeding, however, a convergence factor will be 
selected as the maximum value allowed by Eq. (3.18). Therefore, 
assuming a 1.0 second response window, 

T q = 1.0 
and 


a = 4.605 

Now, in accordance with the procedure outlined in Chapter 3, a 
value of f £ is selected by examining the system transfer function 
shown in Figure 7.1. The transfer function Is effectively zero for 
frequencies greater than 50 Hz, thus 

f c = 50 Hz 

Figure 7.2 shows the exponentially-windowed force spectra at 
frequencies less than f for several Nyquist frequencies. An exami- 

V 

nation of the spectra reveals that as the Nyquist frequency increases, 
the spectra begin to converge to what is assumed to be an exact repre- 
sentation of the exponentially-windowed force spectrum for frequencies 
less than f . For the exponentially-windowed force spectra shown in 
Figure 7.2, the convergence is good for Nyquist frequencies greater 
than 512 Hz. Thus, 
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f H » 512 Hz 

Therefore, having selected T Q and f^ , the renaining 

sampling parameters may be determined using Eqs. (B.l) through (B.4) in 
Appendix B, 

K = 1024 f s = 1024 Hz 

T = 9.7656 * 10" 4 sec f = 1.0 Hz 

o 

The system response may now be determined using Eqs. (3.7) 
through (3.13). The exact solution and the frequency-domain solution 
have been plotted in Figure 7.3a using the sampling parameters deter- 
mined above. The two solutions are nearly identical for f N = 512 Hz . 
However, Figure 7.3b reveals that selecting a Nyqulst frequency which 
is less than that required causes the frequency-domain solution to have 
a significant error. 

B. Selecting a Convergence Factor 

In order to demonstrate the selection of a proper convergence 
factor, the structure used In part A will be subjected to a 1 Hz cosine 
forcing function. Thus, 

p(t) = 987.0 * cos ( ?irt) 

The system response Is desired for t less than 0.6 seconds, thus 

T - 0.6 

Selecting the convergence factor to be the maximum allowed by Eq. 
(3.18) yields 
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a = 7.675 

As in the previous example, the spectrum of the exponentially- 
windowed forcing function shown in Figure 7.4 is examined to determine 
a suitable Nyquist frequency. Using the same criteria as in the 
previous example, the Nyquist frequency is chosen to be 

f = 426.7 Hz 

The remaining sampling parameters may now be determined and the 
response calculated. 

An examination of the exponentially-windowed response, in 
Figure 7.5a, reveals that the convergence factor selected above has 
caused the response to decay to ‘effectively zero' by the end of the 
window. Therefore, the frequency-domain response will not be time- 
aliased. The comparison of the exact and frequency-domain solutions 
shown in Figure 7.5 supports the statement above, as the solutions are 
nearly identical. 

In order to demonstrate the effectiveness of the convergence 
factor in preventing time-domain aliasing, the response will be deter- 
mined using 


a = 1.0 

An examination of the exponentially-windowed response In Figure 
7.6a reveals that the convergence factor has not caused the windowed 
response to be 'effectively zero' by the end of the response window. 
Therefore, as expected, the comparison of the exact and frequency- 









92 



Figure 7.6 (a) Exponentially-windowed response and (b) the 
response with the window removed (a = 1.0). 
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domain solutions in Figure 7.6b is not good, as the frequency-domain 
solution is time-aliased. 

There is, however, a problem with selecting a convergence 
factor which is too large, and which causes the windowed response to be 
'effectively zero' too early in the window. In order to demonstrate 
the problem, let 


a = 30.0 

and let the response window be 1.0 seconds long. Thus, 

T o = L0 ‘ 

The exponentially-windowed response shown in Figure 7.7a 
reveals that the windowed response is 'effectively zero' at t = 0.2 
seconds. Therefore, it is clear that the solution will not be 
time-aliased. However, an examination of Figure 7.7b reveals that the 
frequency-domain solution of the response is not stable near the end of 
the window. Figure 7.7c shows the response for the first 0.8 seconds 
of the window, and it can be seen that the response prior to the region 
of instability compares well with the exact solution. 

Example 2 

This example will demonstrate the substructure coupling proce- 
dures discussed in Chapter 4. The superstructure to be used in this 
example is an 18 degree-of-freedom beam clamped at. each end. 
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Substructure 0 


The beam will be analyzed using the two substructures. 



Substructure 2 


The structure will be loaded at the translational degree of 
freedom of the interface with the following load, 

P(t) 

20 



.01 


t(sec) 
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Assuming the response is desired at only the loaded degree o* 
freedom, the relevant set nay be determined by taking the union of the 
following sets: 

Interface DOF 9,10 

Loaded DOF 9 

Response DOF 9 

Therefore, the relevant set consists of 

Relevant DOF 9,10 

Note that even though the response is not desired at system degree- 
of-freedom 10, it must be determined, since interface degrees of 
freedom are required to be in the relevant set. 

Assuming the first 0.1 second of the response Is desired, 

T q = 0.1 sec 

and the maximum value of the convergence factor is 

a = 46.05 

An examination of a typical transfer function for the system shown in 
Figure 7.8 reveals that the transfer function is not 'effectively 
zero,' even at frequencies greater than 2000 Hz. However, Figure 7.9 
reveals that the exponentially-windowed force spectra is 'effectively 
zero' by 200 Hz, therefore. 


200 Hz 
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An examination of Figure 7.10 reveals that the spectra have nearly 
converged for a Nyqulst frequency of 2560 Hz. Thus, let 

f N = 2560 Hz 

This Nyquist frequency may now be used to obtain the remaining sampling 
parameters and the exponentially-windowed force spectrum. 

An elgensolution may now be performed on each of the substruc- 
tures to obtain the generalized mass, stiffness and modal matrices, 
which are required for the modal formation of the dynamic stiffness 
transfer matrix. The frequency loop and the Inverse transformation 
outlined in Figure 4.3 may now be executed to obtain the time-domain 
response of the desired system degree of freedom. 

The response obtained using the frequency-domain substructuring 
procedure described above is plotted In Figure 7.11 along with an 
accepted time-domain solution. The solutions compare very well, and 
thus the substructuring technique has been verified. 

Example 3 

This example will demonstrate the recovery of superstructure 
natural frequencies from frequency-domain substructure models. The 
superstructure chosen for the example Is the six degree-of-freedom mass 
spring oscillator, 



Substructure 0 

m 3 1.0 k 3 1.0 
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which will be analyzed using the two substructures. 



Substructure 1 



Substructure 2 


The two substructures are joined by one common interface degree of 
freedom. Thus , 


P 


0 

1 


1 


In order to determine the superstructure natural frequencies, 
the fixed Interface natural frequencies of each substructure must be 
obtained. They have been determined and tabulated below. 


Substructure Natural Frequencies 
(rads/sec) 


Substructure 1 

0.4450 

1.2470 

1.8019 


Substructure 2 

0.6180 

1.6180 


The search parameters will be chosen such that the natural 


frequencies will be calculated to within 0.001 rads/sec and the search 
range will be from 0-100 rads/sec. Thus, 
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TOL = 0.C01 (rads/sec) 

“min * 0*000 (rads/sec) 

“mav = 100.0 (rads /sec) 

rod x 

The system natural frequencies may now be obtained using the 
above information and the dynamic stiffnpss transfer matrices of the 
substructures, as directed by the outline in Figure 5.3. 

After an average of 10 polynomial evaluations (i.e. 10 times 
through the bisection loop) per natural frequency, the system natural 
frequencies were determined to be 

Substructure 0 Natural Frequencies 
(rads/sec) 

0.000 

0.518 

1.000 

1.414 

1.733 

1.933 

The characteristic equation of the reduced system Is shown in 

•• 

Figure 7.12. The figure also Illustrates how each of the Sturm se- 
quence Indices changes as a function of frequency. Note that the value 
of S is never greater than one, since pj° is equal to one. From 
the figure, it can be seen that S Is the sum of S and S^ 0 , and 
that the r th natural frequency is located at the frequency where S 


changes from r-1 to r . 
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Example 4 

This final example Is included to demonstrate the applicability 
of Eqs. (6.8) and (6.9) for determining the damped natural frequencies 
of the damped mass-spring oscillator below. 


m 


k k 



c c 


Substructure 0 

1.0 k = 1.0 c a 0.05 


The natural frequencies will be determined using the two substructures. 


k 



c 


Substructure 1 


k 


1 — 

- sAA/— 

I 

Ll 


JU 


c 


Substructure 2 


which have one common Interface degree of freedom, and thus 

p? ■ 1 

The fixed-interface natural frequency of both substructure 1 
and substructure 2 Is (-0.05 + jl.413), since each of the substruc- 
tures are Identical. In order to use the bilinear search procedure 
described in Chapter 6, the range of the imaginary part of the natural 
frequencies has been chosen to be 0-2 rads/sec, and the tolerance as 
0.00001 rads/sec. 
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After an average of 60 polynomial evaluations per root, the 
damped natural frequencies of the superstructure were determined to be 

Superstructure Damped Natural Frequencies 

-.0095 + jO.618 
-.5000 + j 1.413 
-.65*5 + jl.617 

These natural frequencies match precisely the damped natural frequen- 
cies obtained from an exact solution of the equations of motion de- 
scribing the system. 

As noted in Chapter 6, a single example does not explicitly 
prove that Eqs. (6.8) and (6.9) are valid for all systems; however, it 
does demonstrate their applicability in obtaining the damped natural 
frequencies of at least one particular problem. 



Chapter 8 
CONCLUSIONS 

8.1 Solving the Response Problem 

Frequency-domain analysis has been found to be an efficient 
method for finding the frequency response of a structure and has been 
used by structural analysts for many years. The research completed for 
this thesis has also found it to be a suitable method for determining 
the transient response of systems subjected to a wide variety of loads. 
However, since a large number of calculations are performed within the 
discrete frequency loop, the method loses its computational efficiency 
If the load must be represented by a large number of discrete frequen- 
cies. 

It has also been discovered that substructure coupling In the 
frequency domain works particularly well for analyzing structural 
systems with a small number of Interface and loaded degrees of freedom. 
A system with these restrictions will have relatively small matrices 
within the frequency loop, which increases the efficiency of the 
frequency-domain procedure. 

The ability to describe large complex substructures by the 
dynamic stiffness and loads at only the interface degrees of freedom 
makes substructure coupling In the frequency domain very attractive for 
structural modification problems as well. 
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8.2 Determining the Natural Frequencies 

In an attempt to determine the generality of frequency-domain 
analysis, it was discovered that substructure coupling in the frequency 
domain can lead to an efficient method of obtaining natural frequencies 
of undamped structures. The efficiency, however, decreases as the 
number of interface degrees of freedom increases. 

It has also been found that the damped natural frequencies of a 
system may be determined using frequency-domain techniques. However, 
the algorithm used to determine these natural frequencies has not been 
proven explicitly. 

8.3 Directions for Future Research 

Over the past two decades substructure coupling In the tine 
domain has played an important role In the analysis and design of 
structures. In order for substructure coupling in the frequency domain 
to play a similar role, it will be necessary for the following items to 
be investigated: 


- A more In-depth study on the effects of using 
a convergence factor In transient response 
analyses. 

- A solution procedure for analyzing response 
to random excitation. 

- The effects which modal truncation at the 
substructure level has on the system response 
and natural frequencies. 



A proof of Eqs. (6.8) and (6.9), which were 
used to determine the damped natural fre- 
quencies of the superstructure. 
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With these topics clarified, substructure coupling in the 
frequency domain could play an important role in the design and analy- 
sis of future structural systems. 



APPENDIX A. TRANSFORMING THE EQUATIONS OF MOTION 


A.l Unilateral Fourier Integral Transform 

The unilateral Fourier transform pair chosen for this thesis 
may be written as 


YU) = f y(t) e’ jwt dt (A.l) 

Jo 

y(t) - ^ f YU) e jwt da) (A. 2) 

In order for the unilateral Fourier Integral transform of y(t) to 
exist, and for y(t) to be recoverable from Its transform, the 
following conditions must be satisfied: 

a) In every finite Interval, y(t) must be bounded and have a 
finite number of maxima and minima, and a finite number of 
discontinuities. 

b) The integral 

I!! L |y(t)| dt 5 f |y(t) | dt (A.3) 

converges (l.e. the Integral approaches some finite limit 
as A approaches infinity). 

Since both the system excitation and response will be 
transformed to the frequency domain, each must satisfy the conditions 
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stated above. The first condition will be satisfied for any stable 
system subjected to a physically-realizable excitation. The second 
condition may be satisfied by many functions. However, this discussion 
will assume the excitation and response are identically zero beyond 
some finite value of t . This assumption will ensure convergence of 
the Integral given by Eq. (A. 3) (Ref. LePage). 

Therefore, with the transformability conditions satisfied, the 
system equation of motion, Eq. (2.1), may now be transformed to the 
frequency domain using Eq. (A.l). 


J~ [[ml{x} + [c](x> + [k]{x>] e"^ ut dt 

* f (f(t)} e' jwt dt 
In 


(A. 4) 


The Integration of the first two terms by parts yields 


[m]{x> 

e’ ja>t dt = [m]{x> 

e -j«t 

r 




0 

+ 

,1to[ml{x> e“^ ut (q - 

to 2 [m] 

j (x> e" jut dt 

j 0 

[c]{x> 

e’ j “t dt = [c]{x> 

e -J“t 

r 




0 

+ 

Me] f tx> 5 3 “* 

dt 



(A.5) 


(A. 6) 


and combining terns yields 
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[-<o 2 [m] + j«[c] + [k! ] (x) dt 


= f {f(t)> e‘ jut dt - [m]{x> e" jtot 1" 
J 0 0 

(A.7) 

- [Mm] + [c]](x> e- jut |“ 


The evaluation of the last two terms on the right hand side of 
Eq. (A.7) at the lower limit Is accomplished by using the Identity 

e -M , = 0 = , 

e 't a 0 e 1 

(A. 8) 

Although e”^ Is not defined at t = » , Itls known that 
| e -J“"| s i # and therefore the terms may be evaluated at their upper 
limit by recalling the assumption of finite-duration response. That is 
to say. 

(xH t _ - (0} 

(A. 9) 

<x)l t „ - ( 0 ) 

(A. 10) 

Therefore, Eq. (A.7) may be simplified and written 


[-<u 2 [m] + jw[c] + [k]] (X(w)> = (F(u)> + [m](x 0 > 

(A. 11) 


+ [Mm] + [c]] (x 0 ) 


where 


FU) - 

1 (f(t)> dt 

Jo 

(A. 13) 

<v ■ 

<x(t=0)} 

(A. 14) 

<*„> = 

{x(t=0)> 

(A. 15) 


Eq. (A. 11) represents the system equations of motion In the frequency 
domain, and is subject to the restriction that the excitation be 
nonperiodic and of finite duration. It is also restricted to those 
systems which are damped, in order that the response will also be of 
finite duration. 

A. 2 Unilateral Fourier Integral Transform with Convergence 
Functions 

The unilateral Fourier integral transform pair with convergence 


functions may be written as 


YU) 

* f* y(t) e"^ wt dt 

Jo 

(A. 16) 

y(t) 

• f YU) e jwt do, 

(A. 17) 


where 
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In order for the unilateral Fourier integral transform of v(t) to 
exist, and for y(t) to be recoverable from its transform, the 
following conditions must be satisfied. 

A 

a) In every finite interval, y(t) must be bounded and have a 
finite number of maxima and minima, and a finite number of 
discontinuities. 

b) The integral 

f |y(t)|e" at dt (A. 19) 

J 0 

converges . 

As in Appendix A.l, both the system excitation and response 
must satisfy the conditions stated above, and the first condition is 
satisfied for any stable system subjected to a physlcally-realizable 
excitation. However, the second condition will not restrict the class 
of problems which can be worked, so long as a value of the convergence 
factor, a , may be found which is greater than zero and causes Eq. 
(A. 19) to converge. 

The equation of motion, Eq. (2.1), may now be transformed using 
Eq. (A. 16). Thus, 


[[m]{x> + [c](x> + [k]{x>] 


e" at e"** dt 


= f {f(t)> e’ at e~i ut dt 

Jo 


(A. 20) 


and simplifying yields 
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J~ [[m]{x> + Cc](x> + [k]{x>] e" st dt 


f (f(t)} e~ st dt (A. 21) 

> 0 


where the complex frequency variable, s , is given by 


s = a + ju 


(A. 22) 


The integration by parts of the first two terms yields 


[[m]s 2 + [c]s + [k]] {x> e“ st dt 

|* (f(t)} e _st dt - [m]{x> e _st |" - [ [m]s + Tc]] (x> e 


(A. 23) 


-St |® 
0 


Now, the last two terms on the right hand side of Eq. (A. 23) may be 
evaluated at the lower limit using the identity 

e’ st | t=0 ■ e° = 1 (A. 24) 


and, at the upper limit, since a > 0 



Therefore, Eq. (A. 23) may be simplified and written 
[[m] s 2 + [c] s + [k]]{X(s)> 

= <F(s) > + [m] {x Q } + [[m] s + [cl] {x Q } 


(A. 25) 


(A. 26) 


where 
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{X ( s ) > = f (x) e“ st dt 

J 0 

(A. 27) 

(F(s ) > = f {f(t)> e" st dt 
j 0 

(A. 28) 

(x Q } = {x (t=0)} 

(A. 29) 

(x 0 > = (x (t=0)} 

(A. 30) 


Equation (A. 26) represents the system equations of motion in the 
complex frequency domain. The only restriction on the system, is that 
a convergence factor exists such that Eq. (A. 19) converges. 

A. 3 Complex Fourier Series 

Assume that the system described by Eq. (2.1) is subjected to a 
periodic excitation which may be represented by 

{f(t)> = {F} e jnt (A. 31) 

Now, assume the response will be represented In a similar form 

{x(t)> = (X> e jnt (A. 32) 

Equation (A. 31), along with Eq. (A. 32) and Its derivatives may now be 
substituted into Eq. (2.1) to yield 

jilt _ gjilt 


[-fi 2 [m] + jn [c] + [k]J (X) e 


(A. 33) 
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or in a simplified form as 

[G(fl)] (X) = (F> (A. 34) 

Equation (A. 34) may now be solved for (X) , which is then substituted 
into Eq. (A. 32) to obtain the steady-state response of the system to 
the periodic excitation (f(t)} 


(x(t)> = [G(n)]“ 1 (F> e jnt (A. 35) 


For physical ly-reallzable periodic excitations which are 
composed of more than one frequency, the excitation can be represented 
by the complex Fourier series transform pair 

{f(t)J = Z (Ffo). )) e j<Ukt (A. 36) 
'l k=-» K 


fT+T- -ju.t 

fF(co k )> = 1 {f(t)J e k dt 


(A. 37) 


where 



(A. 38) 


and Tj is the fundamental period of the excitation. 

For linear systems, the principal of superposition may be 
invoked to obtain the total response 



£ 

k=-» 


[G(u k )] _1 {F(o> k )> e K 


(x(t ) > 


(A. 39) 



where 
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[G(w k )l = -o> k 2 [m] + jto k [cl + [k] (A. 40) 


In order that a frequency-domain equation may be written In a 
form similar to those In Appendices A.l and A. 2, Eq. (A. 39) will be 
transformed to the frequency domain and written 


[-*o k 2 [m] + J« k [c] + [k]] {X(a> k )} = {F(u k )> (A.41) 


k s _CO 0 

I • • • y 


where 


ft+T. -ju.t 

tX((u k )} = J 1 (x(t)> e K dt (A. 42) 

and {F(co k )} Is given In Eq. (A. 37). 

Therefore, Eq. (A.41) represents the equation of motion In the 
frequency domain which describes the steady-state response of a system 
subjected to periodic excitation. 



APPENDIX B. DISCRETE FOURIER TRANSFORM 


In the discussion of the Discrete Fourier Transform (DFT) and 
the Fast Fourier Transform (FFT) in Chapter 3, the following notation, 
parameter relationships and terminology are often utilized. 

Notation 


f s - 

sampling frequency 


f N ' 

Nyqulst frequency 


f o - 

frequency resolution 


T „ - 

length of the time record 


T - 

time resolution 


K - 

number of discrete samples 


Relationships 



f N ■ 

1 

TT 

(B.l) 

f s ' 

2f N 

(B.2) 

T o ■ 

1 

(B.3) 


T 


T - 

0 

IT 

(B.4) 


Terminology 

Frequency-Domain Aliasing - Frequency-domain aliasing occurs when the 
sampling frequency is less than twice the highest frequency contained 
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in the continuous function being transformed. It is characterized by 
the energy of the frequencies which are greater than f^ , folding 
back into the lower frequencies. This increased energy in the lower 
frequencies of the force spectrum will cause the response spectrum to 
also have higher energies in the low frequencies. 

Time-Domain Aliasing - Time-domain aliasing occurs when the frequency 
resolution is not fine enough to allow the Inverse transform of a 
frequency function to be effectively zero by the end of the time record 
T . It is characterized by the nonzero part of the function extending 
beyond T Q , adding back into the function at the beginning of the 
time record. 

Leakage - Leakage occurs when the time record does not contain an 
Integer number of fundamental periods of a periodic function. It is 
characterized by energy from the fundamental frequency showing up in 
frequencies adjacent to the fundamental frequency. 
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